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Abstract 


Numerical  methods  are  used  to  study  the  growth  of  waves  of 
finite  amplitude  on  a  pair  of  parallel  infinite  vortices.  The  vortices 
are  treated  as  lines  except  in  so  far  as  the  detailed  structure  of  the 
core  is  needed  to  remove  consistently  the  singularity  in  the  line 
integrals  for  the  velocities  of  the  vortices.  It  is  shown  that  the 
vortices  eventually  touch  and  the  shape  of  the  wave  at  this  instant 
is  calculated.  The  wave  is  quite  distorted  at  this  instant,  but  it  is 
shown  that  its  gross  properties  are  given  roughly  by  linear  theory. 


The  trailing  vortex  system  behind  a  large  aircraft  contains  a 
considerable  amount  of  energy  and  presents  a  hazard  to  other,  particularly 
smaller,  aircraft.  The  recognition  of  this  problem  has  lead  to  attempts, 
both  theoretical  and  experimental,  to  estimate  the  life-time  of  the 
trailing  vortex  system. 

To  do  this  theoretically  one  must  isolate  the  mechanism  which  is 
responsible  for  destroying  the  vortex  trail.  Turbulent  diffusion  of  the 
vorticity  is  too  gradual  to  matter  and  only  two  processes  seem  able  to 
cause  rapid  disintegration  of  the  vortex  tystem,  vortex  bursting  (Bisgood 
et  al  1971)  and  the  growth  of  waves  of  large  amplitude. t  This  latter 
effect  can  often  be  observed  when  the  vortex  system  of  a  high-flying 
aircraft  is  rendered  visible  by  condensation.  Wavy  disturbances  grow 
on  both  trailing  vortices  and  rapidly  reach  an  amplitude  such  that  the 
two  vortices  touch  at  the  points  of  maximum  inward  displacement.  The 
trail  then  breaks  up  into  a  sequence  of  distinct  rings  and  soon  ceases  to 
be  visible.  Several  systematic  observational  studies  of  this  "looping” 
process  have  recently  been  made  (Rose  and  Dee  1963,  Crov'  and  Murman 
1970,  Olsen  1971,  Widnall,  Bliss  and  Zalay  1971). 

It  is  worth  emphasizing  that  the  mechanism  by  which  the  original 
vortices  break  and  rejoin  to  form  rings  is  not  understood.  Nor  is  it 
clear  why  the  vortex  rings,  which  are  stable  on  classical  theory,  should 
disappear  so  rapidly  (if  indeed  they  do,  since  care  must  be  taken  in 

^At  high  altitudes  \arge  stable  density  gradients  can  occur,  and  Scorer 
and  Davenport  (1970)  have  suggested  that  buoyancy  effects  can  then  lead 
to  rapid  disintegration. 
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interpreting  observations  which  depend  on  the  retention  by  the  vortex 
of  smoke  particles  or  water  droplets).  However,  the  present  paper  is 
concerned  only  with  the  motion  of  the  vortices  up  to  die  instant  at  which 
they  touch. 

Theoretical  work  on  the  growth  of  waves  was  initiated  by  Crow 
(1970).  Crow  represented  the  trailing  vortex  system  as  a  pair  of  infinite 
parallel  vortex  filaments  of  circulation  +.  r  having  circular  cross  sections 
of  radii  a  small  compared  to  their  separation  b.  Small  perturbations  in 
the  form  of  plane  waves  were  shown  to  be  unstable,  provided  only  that  the 
waves  were  sufficiently  long,  the  time  of  growth  being  of  order  b2/r. 

In  the  mode  of  most  rapid  growth  the  disturbances  lie  in  planes  inclined 
at  about  45*  to  the  vertical  which  are  oriented  so  as  to  meet  below  the 
undisturbed  positions  of  the  vortices  and  the  waves  themselves  are 
symmetric  in  the  sense  that  the  crests  of  the  wave  on  the  left-hand  vortex 
are  directly  opposite  the  crests  of  the  waves  on  the  right-hand  vortex. 

The  wavelength  of  this  mode  of  maximum  growth  rate  depends  weakly 
on  a/b  ancl  for  the  value  taken  from  the  simplified  theory  of  Spreiter 
and  Sacks  (1951),  Crow  showed  that  \/b  =  8.5. 

According  to  linear  theory  the  troughs  of  the  waves  will  descend 
in  their  respective  planes  until  they  meet,  so  that  an  explanation  of  the 
first  stage  of  the  looping  process  ie  to  hand.  However  the  waves  are  of 
amplitude  comparable  to  their  length  when  they  touch,  so  that  it  is  clearly 
desirable  to  examine  non-linear  effects. 

A  first  attempt  has  been  reported  by  Hackett  and  Theisen  (1971) 
and  it  was  shown  that  while  the  wave  is  modified  by  non-linear  effects, 
the  troughs  do  eventually  meet  in  the  manner  suggested  by  linear  theory. 
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However  their  numerical  method  did  not  permit  the  effect  of  axial  flow 
in  the  vortex  filaments  or  differing  filament  radii  to  be  studied  and  it 
was  thought  worthwhile  to  present  a  numerical  method  in  which  these 
effects  could  be  allowed  for. 

The  assumptions  made  are  precisely  those  of  Crow's  theory.  The 
fluid  is  regarded  as  inviscid  and  incompressible  and  has  uniform  density, 
while  the  vortex  system  is  modelled  as  described  above.  To  simplify  the 
calculation  still  further,  the  vortices  are  treated  as  being  of  zero  cross 
section  and  the  velocity  field  due  to  either  vortex  is  calculated  from  the 
Biot-Savart  line  integral.  The  neglect  of  viscous  effects  means  that 
these  vortex  lines  move  with  the  fluid,  and  this  leads  to  a  pair  of  coupled 
integro-differential  equations  for  the  equations  of  the  vortex  lines. 

A  difficulty  occurs  because  the  Biot-Savart  integrals  diverge  on 
the  vortex  lines  themselves,  but  this  can  be  dealt  with  by  means  of  a 
"cut-off.  "  The  choice  of  this  cut-off  depends  on  the  detailed  structure 
of  the  vortex  filament,  particularly  its  radius  and  the  radial  distribution 
of  vorticity  and  axial  flow.  It  is  only  through  the  cut-off  that  these 
details  influence  the  calculation. 

These  matters  are  described  more  fully  in  §2  and  in  §3  some 
transformations  required  before  efficient  numerical  integration  is  possible 
are  described. 

In  §4  tne  results  of  numerical  integration  are  described  for  the 
case  already  discussed  on  linear  theory  by  Crow.  There  is  no  axial 
flow  and  the  vorticity  is  uniform  inside  the  filament,  which  has  the 
radius  given  by  the  Spreiter  and  Sacks  theory.  Initial  plane  waves  of 
small  amplitude  and  wavelength  and  orientation  corresponding  to  Crow's 
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most  unstable  mode  are  prescribed.  The  subsequent  growth  o£  the 
disturbance  is  then  followed  by  step-wise  integration  of  the  integro- 
differential  equations. 

The  wave  remains  almost  plane,  bat  is  somewhat  distorted  from 
a  sinusoidal  form  by  the  time  the  troughs  meet. 

In  view  of  the  lack  of  reliable  information  about  the  distribution 
of  vorticity  and  axial  velocity  in  the  filament  in  the  real  case  it  was  not 
thought  to  be  worthwhile  to  carry  out  further  calculations  at  present, 
though  the  program  can  deal  with  the  general  case. 
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$2.  The  equations  of  motion 

Rectangular  axes  moving  with  the  undisturbed  vortices  are  employed, 
Oz  being  parallel  to  the  vortices  and  midway  between  them  and  Ox  being 
perpendicular  to  the  plane  defined  by  the  undisturbed  positions  of  the 
vortices.  For  an  aircraft  flying  horizontally  Oz  will  be  in  the  direction 
of  flight  and  Ox  will  be  vertically  upwards. 

The  vortices  will  be  represented  parametrically  as  x^(v,t)  and 
x^(v,t)  where  v  is  a  parameter  chosen  so  that  vs  const,  represents 
always  the  same  fluid  particle.  Then,  bearing  in  mind  that  the  fluid 
at  large  distances  has  a  velocity  i  relative  to  the  moving  axes,  the 

equation  of  motion  of  the  left-hand  vortex  is 


/0XT  v 

.  X 

dv 

Ut  ) 

4ir 

T.  A  ... .  ,.ij 

0V  |xl(v0)-xl(v)|3 


.  -X 

4ir 


oo  3xn 

J  dv  A 

-oo 


feL(v°)^R(v))dv 

|£l(v°)"2r(v)|J 


(2.1) 


and  that  of  the  right-hand  vortex  is 
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The  notation  £  implies  that  some  artifice  is  used  to  make  the  integrals 
representing  the  self  induction  finite  at  v=v0  or  v=v2. 
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Two  methods  of  dealing  with  the  divergence  of  the  integral  for 
the  self-induced  velocity  have  been  employed.  The  simplest  method 
and  that  used  by  Crow  himself  is  to  remove  an  interval  (v0-  e,  v0  -c ) 
from  the  range  of  integration,  where  the  small  quantity  e  is  chosen 
so  that  the  arc  length  of  the  portion  of  the  vortex  corresponding  to 
v0-e<  v  <  v0  +  c  bears  a  fixed  ratio  26c  to  the  core  radius  a(v0)  at 
the  point  in  question.  The  quantity  6c  is  chosen  so  that  the  velocity 
given  by  the  "cut-off"  integral  agrees  with  that  given  by  the  exact  solution 
for  a  circular  vortex.  Remarkably  this  choice  of  can  be  shown  to 

give  the  correct  self-induced  velocity  for  a  vortex  of  any  shape  provided 
that  terms  of  0(  — )  are  neglected,  where  p  is  the  radius  of  curvature  of 
the  vortex  (Widnall,  Bliss  and  Zalay  1971,  Moore  and  Saffman  1971). 

This  method  of  cut-off  is  not  flexible  enough  for  numerical  work, 
because  it  would  be  awkward  to  remove  an  interval  from  the  line  integral 
which  was  not  terminated  by  grid^points  of  the  finite  difference  scheme. 

It  may  be  remarked  at  this  point  that  dealing  with  the  infinity  by  simply 
omitting  the  grid  point  at  which  the  integral  is  infinite,  as  was  done  by 
Hackett  and  Theisen  (1971),  is  equivalent  to  a  choice  of  core  radius 
dictated  by  the  mesh-size  of  the  finite  difference  scheme. 


A  suitable  method  of  cut-off  is  that  suggested  by  Rosenhead  (1930), 


and  consists  of  replacing  the  denominator  in  the  integral  for  the  self 

Vz 


where  the  small  length 


induced  velocity  in  (2.1)  by  J  (?Sijvo)"JSjv))2  +  H2 
p.  bears  a  fixed  ratio  26R  to  the  core  radius  a(v0)  and  similarly  for  (2.  2). 
Rosenhead  merely  remarked  that  p  is  0(a)  but  by  comparing  the  solution 
using  a  cut-off  integral  with  the  exact  solution  for  a  ring,  as  generalized 
by  Widnall  et  al  (1971),  it  can  be  shown  that 


2.3 


1o«26r  =  -i  -  rTjT)  /  r»(.)  *  +  /‘  »w»d.  .  (2.3) 

0  0 

where  r(s)  is  the  circulation  and  W(s)  the  axial  velocity  at  a  distance 
s  from  the  centre  of  the  vortex  filament.  For  W  =  0  and  r(s)  =  s*/al  * 
so  that  there  is  no  axial  flow  and  uniform  vorticity  inside  the  filament 
(2.  3)  reduces  to 


log  26r  =  -  |  (2.4) 

For  the  reasons  described  in  §1  this  is  the  only  case  considered  in  detail. 
However  one  is  now  faced  by  a  more  difficult  problem,  which  is  how  the 
core  radius  itself  is  to  be  determined,  since  until  a(v0)  is  known  p  cannot 
be  determined.  The  vortex  filaments  will  stretch  as  the  waves  grow  so 
that,  since  their  volume  is  conserved,  their  cross-sectional  area  will 
certainly  have  to  change. 

This  question  has  been  studied  elsewhere  (Moore  and  Saffman 

1971)  and  it  has  been  shown  that  a  is  sensibly  uniform  along  the  vortex, 

being  just  that  function  of  time  only  which  conserves  the  volume  of  the 

filament.  Essentially  this  is  due  to  the  fact  that  changes  of  cross - 

p 

sectional  area  propagate  along  the  vortex  with  speed  O(-)  ,  whereas 
changes  of  curvature  or  torsion  travel  with  speed  C -  log 
Thus  variations  of  cross-sectional  area  are  smooth  d  <<  t  in  a  time 
short  compared  to  the  times  of  interest  here. 

This  is  a  great  simplification  and  means  that  the  radius  of 
the  filament  is  easily  determined  from  the  instantaneous  length  of  the 
filament. 
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It  is  now  a  simple  matter  in  principle  to  determine  the  evalution  of  any 
initial  disturbance  g^(v, 0)  ,  (v,  0)  to  the  vortex  pair,  by  integrating 

(2.1)  and  (2.2)  forward  in  time,  the  length  of  the  filaments  being 
calculated  at  each  time  step  in  order  to  obtain  the  appropriate  value 
of  p.  However  some  transformations  of  (2.1)  and  (2.2)  which  take 
account  of  the  periodic  nature  of  the  disturbance  greatly  decrease  the 
amount  of  computing  necessary  and  these  are  described  in  the  next 
section. 


Mliiii’iinrrif  *t-  - 
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§3.  Reduction  to  a  form  suitable  for  numerical  integration 

Only  symmetric  disturbances  are  to  be  considered  so  that 


xL(v,t)  =  Xj^V.t)  , 
yL(v,t)  =  -yR(v,t)  , 
*L(v,t)  =  zR(V,t)  . 

J 


(3.1) 


Hence  only  (2.1)  need  be  integrated  and  (2.2)  is  replaced  by  the  symmetry 
conditions  embodied  in  (3.1). 

Next  advantage  is  taken  of  the  fact  the  disturbances  are  periodic, 
which  means  that  the  parameter  v  can  be  chosen  so  that  for  any  integer  n 


XL(v  +  n,  t)  =  k  n  X  +  xL(v,t)  . 


(3.2) 


Here,  X  is  the  wavelength  of  the  disturbance,  which  does  not  change  as 
the  wave  grows,  and  v  has  been  chosen  so  that  0  <  v  <  1  corresponds 
to  one  wavelength. 

Thus  one  need  only  follow,  say,  the  portion  -j  ^  v  <  \  to 
obtain  Xjjv.t)  and  hence,  by  (3.1)  ,  .  But  the  integrals  in  (2.1) 

converge  rather  slowly  and  have  large  integrands,  so  that  further  steps 
are  needed  to  make  accurate  numerical  integration  feasible. 

The  first  integral  in  (2.1)  can  be  decomposed  into  the  form 


8£l  (2L(Vo)'2L(v))dv  v  x+l  8xL  feL(v®)'^L(v)'2m^V)dv 

y  17  a- - “T" +  L  y  IT" a  , - 77  0.3) 

■'  I*L(v°)_2l<v)I  m=  _I  I^L(Vo)'XL(v)_2m~Xl 
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'where  the  second  term  arises  from  summing  over  successive  portions  of  the 
dsturbance  of  length  2\.  If  v0  is  restricted  to  lie  in  the  interval  {-£  ,  £)  , 
all  the  integrands  in  the  second  term  are  bounded  t  and  no  cut-off  is 
required.  Thus  the  integral  can  be  written  as 


a  £<v»vo>dv 


-i 


8v 


|A<  'j  *  I*'  TT  A  (A(v.v0)P(A)-kQ(A))dv 
f£(v.v0)|  -> 


(3.4) 


where 


A(v,v0)  =  xL(v0)  -  xL(v)  (3.5) 

and 

P(u) 

Q(u) 

The  sums  of  the  infinite  series  in  (3.6)  do  not  seem  to  be  expressible  in  terms 
of  elementary  functions,  but  they  can  be  computed  numerically.  In 
practice  a  table  of  values  over  a  discrete  set  of  uz  and  u  •  k  was 

»V 

computed  and  stored  before  time  integrations  were  started,  the  values 
of  P  and  Q  being  calculated  for  the  actual  u  arising  in  the  time  integration 
by  bi-variate  interpolation. 

^This  is  why  the  summation  is  over  two  waves  at  a  time.  I  am 
indebted  to  Professor  P.  G.  Saffman  for  this  step. 
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(3.6) 
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This  use  of  the  periodicity,  which  was  suggested  by  Rosenheads 
(1930)  treatment  of  the  two-dimensional  vortex  sheet,  means  that  the 
range  of  integration  is  reduced  from  (-ca,oo)  to  (-1,1)  so  that  high 
accuracy  can  be  achieved  in  the  numerical  integration  with  a  tractable 
number  of  grid  points. 

The  second  integral  in  (2.1)  can  be  similarly  transformed  and 
the  same  functions  P  and  Q  arise. 

However  the  integral  in  the  first  term  in  (3.4),  though  every¬ 
where  finite,  is  large  near  v0  and  this  would  cause  loss  of  accuracy  in 
evaluating  the  integral.  It  is  easy  to  show  that,  when  v  ~  v0 , 


so  that  if  the  first  term  in  (3.  4)  is  written  in  the  form 


then  the  integrand  in  the  first  integral  is  0(1)  everywhere,  while  the  second 
integral  is  elementary. 

Finally  it  may  be  remarked  that,  in  view  of  the  uniformity  of 
the  cross  section  and  the  conservation  of  volume, 
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i*  *  26R»o{i  fv‘  .  t».») 

-'A 

where  a0  is  the  uniform  radius  of  the  vortex  filaments  in  the  undisturbed 
state. 

It  is  now  a  straightforward  matter  to  follow  numerically  the 
evolution  of  a  prescribed  initial  disturbance.  The  interval  (-1.1)  was 
divided  into  2N  portions  by  2N+1  equally  spaced  grid  points.  The  spatial 
derivatives  were  calculated  using  four-point  centered  differences  and 
Simpson's  rule  was  used  to  evaluate  the  integrals. 

The  integration  forward  in  time  was  effected  by  the  fourth-order 
Runge-Kutta  formula,  this  method  being  used  because  of  its  stability. 

Clearly  dimensionless  variables  are  preferable  in  the  numerical 
work  and  variables 

x'  =  x./X 

^  fx  I  4 

and  (3.10) 

t'  =  t(r/4ff\2) 


were  found  to  be  most  convenient. 


The  results  of  a  typical  calculation  are  described  in  §4. 


§4.  Numerical  Results 

Crow  (1970)  has  shown  that  the  natural  time  scale  for  the  growth 
of  waves  is  2irb*/r  ,  which  is  just  the  time  taken  for  the  undisturbed 


vortices  to  descend  a  distance  equal  to  their  separation.  Thus  to  present 

$ 

the  results  a  dimensionless  time  t  defined  by 


(4.1) 


will  be  employed. 

The  trailing  vortices  are  assumed  to  have  the  radius  given  by 
Spreiter  and  Sacks  (1951)  so  that 


a/b  =  .098 


(4.2) 


It  must  be  stressed  that  it  is  far  from  certain  that  this  is  realistic  for 
real  trailing  vortices  (see,  for  example  McCormick,  Tangier  and  Sherrieb 
1968)  but  as  no  firm  value  has  emerged  from  experiment  it  seemed  best 
to  adhere  to  (4.2).  In  any  case  the  results  are  not  very  sensitively 
dependent  on  a/b. 

Crow's  linear  theory  shows  that  for  this  value  of  a/b  ,  the  most 

o 

unstable  mode  has  wavelength  8.5b  and  lies  in  a  plane  inclined  at  47.5 
to  the  horizontal.  It  was  decided  to  study  the  evolution  of  this  mode 
from  an  initial  semi-amplitude  of  0.05b. 

* 

Trial  and  error  showed  that  a  time  step  At  =0.1  gave  adequate 

accuracy  and  the  value  N  =  40  (that  is  40  points  per  wave)  was  used. 

$ 

However  at  t  =2.4  numerical  instability  sets  in  at  the  troughs  of  the  waves, 


4.2 


presumably  because  the  distance  between  the  grid  points  was  not  small 
compared  to  the  separation  of  the  vortices  and  because  the  curvature 
of  the  wave  is  largest  here.  To  cure  this  and  to  cater  for  the  more 
rapid  changes  consequent  on  the  proximity  of  the  troughs  the  calculation 
was  stopped  at  t  :  2,2  and  restarted  with  N  =  80  and  At  =  0.25.  This 
removed  the  instability. 

if 

At  t  a  2.475  the  troughs  of  the  vortices  were  only  .1876  apart, 
which  implies  that  the  two  vortices  are  touching.  Clearly  the  results  of 
a  calculation  which  treats  the  separation  of  the  vortices  a&  large  compared 
to  their  radii  must  be  viewed  with  skepticism  at  this  stage.  However,  a 
separate  numerical  study  of  the  motion  of  cylindrical  vortices  in  two 
dimensions  in  which  core  size  was  allowed  for  showed  that  the  Biot-Savart 
formula  gives  roughly  the  correct  velocity  even  when  the  cores  are 
touching.  Thus,  while  the  cores  are  very  distorted,  which  would  affect 
the  internal  structure  and  thus  change  the  cut-off,  it  is  possible  that  the 
approximations  on  which  the  present  paper  is  based  are  adequate  even 
when  the  cores  are  close. 


In  any  case,  there  is  no  point  in  integrating  beyond  the  approximate 

j*C 

instant  of  touching  so  the  calculation  was  terminated  at  t  s  2.475. 

Figure  1  shows  log  10  [” B (t^j/B (0)1  as  a  function  of  t*  ,  where 


B 


y  -  y 

7  max  7  min 

y  +  y 

7  max  7  nun 


(4.2) 


is  a  convenient  measure  *  (introduced  by  Crow  and  Murman  1970)  of  the 


*This  quantity  can  be  directly  measured  from  photographs  taken  from 
vertically  below  the  trail. 


4.3 


growth  of  the  waves.  According  to  Crow's  linear  theory 

log  jo  ^B(t*)/B(0)j  =  .  368 1*  (4.  3) 

and  this  straight  line  is  shown  for  comparison.  The  agreement  is  excellent 
for  small  times  and  indeed  linear  theory  is  adequate  over  nearly  the  whole 
range  of  values  of  t  . 

Figure  2  shows  a  plan  view  of  the  centre-lines  of  the  vortices  at 

i!jH 

the  approximate  instant  of  touching  t  =  2.475,  while  Figure  3  shows  a  side 
elevation  at  that  time.  For  the  reasons  given,  it  is  hard  to  assess  the 
relevance  to  the  actual  situation  of  these  results,  but  there  is  some 
similarity  to  observed  large  amplitude  waves  in  that  the  curvature  is  longer 
at  the  trough  than  at  the  crest.  The  development  of  large  curvature  in  the 
trough  suggests  that  even  more  grid  points  would  have  been  desirable 
the  final  stages  of  the  calculation,  but  in  view  of  the  dubious  relevance  of 
this  part  of  the  calculation  such  a  refinement  was  not  attempted. 

The  end  elevation  is  shown  in  Figure  4  and  reveals  that  the  wave 
remains  practically  plane,  though  the  plane  is  not  quite  the  original  one. 

The  increase  of  length  of  the  wave  was  found  to  be  only  5%  at  t  *  2.475. 

The  vortex  stretches  most  at  the  trough  where  the  distance  between  two 
fluid  particles  initially  a  distance  6z  apart  increases  to  1.6  6z  at  the  instant 
of  touching.  The  vortex  contracts  at  the  crest  where  the  distance  between 
two  fluid  particles  initially  a  distance  5z  about  decreases  to  0.  9  6z.  The 
differences  reflect  the  marked  asymmetry  between  trough  and  crest. 
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Figure  Captione 


Figure  1.  The  function  log10  ^B(t*)/B(0)J  as  a  function  of  t*  obtained  from 
linear  theory  and  from  the  numerical  solution  (circled  points). 

Figure  2.  The  plan  view  of  the  wave  at  t  =  2.475,  the  approximate  instant 

of  touching.  The  initial  core  diameter  is  shown  on  the  same  scale. 

«  >5* 

Figure  3.  The  side  elevation  of  the  wave  at  t  s  2.475. 
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Figure  4.  The  end  elevation  of  the  wave  at  t  =  2.475. 


Figure  2.  The  plan  view  of  the  wave  at  t*  =  2.475,  the  approximate  instant 
of  touching.  The  initial  core  diameter  is  shown  on  the  same  scale. 


Figure  4.  The  end  elevation  of  the  wave  at  t*«  2.475. 


